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There are a number of situations in which several signals are simultaneously recorded in complex 
systems, which exhibit long-term power-law cross-correlations. The multifractal detrended cross- 
correlation analysis (MF-DCCA) approaches can be used to quantify such cross-correlations, such 
as the MF-DCCA based on detrended fluctuation analysis (MF-X-DFA) method. We develop in 
this work a class of MF-DCCA algorithms based on the detrending moving average analysis, called 
MF-X-DMA. The performances of the MF-X-DMA algorithms are compared with the MF-X-DFA 
method by extensive numerical experiments on pairs of time series generated from bivariate fractional 
Brownian motions, two-component autoregressive fractionally integrated moving average processes 
and binomial measures, which have theoretical expressions of the multifractal nature. In all cases, 
the scaling exponents h^y extracted from the MF-X-DMA and MF-X-DFA algorithms are very close 
to the theoretical values. For bivariate fractional Brownian motions, the scaling exponent of the 
cross-correlation is independent of the cross-correlation coefficient between two time series and the 
MF-X-DFA and centered MF-X-DMA algorithms have comparative performance, which outperform 
the forward and backward MF-X-DMA algorithms. For two-component autoregressive fractionally 
integrated moving average processes, we also find that the MF-X-DFA and centered MF-X-DMA 
algorithms have comparative performance while the forward and backward MF-X-DMA algorithms 
perform slightly worse. For binomial measures, the forward MF-X-DMA algorithm exhibits the best 
performance, the centered MF-X-DMA algorithms performs worst, and the backward MF-X-DMA 
algorithm outperforms the MF-X-DFA algorithm when the moment order g < and underperforms 
when q > 0. We apply these algorithms to the return time series of two stock market indexes and 
to their volatilities. For the returns, the centered MF-X-DMA algorithm gives the best estimates 
of hxy{q) since its hxy{2) is closest to 0.5 as expected, and the MF-X-DFA algorithm has the 
second best performance. For the volatilities, the forward and backward MF-X-DMA algorithms 
give similar results, while the centered MF-X-DMA and the MF-X-DFA algorithms fails to extract 
rational multifractal nature. 

PACS numbers: 05.45.Tp, 05.45.Df, 05.40.-a, 89. 75. Da, 89.65.Gh 



I. INTRODUCTION 

Natural and socio-economic systems are usually com- 
plex systems from which macroscopic statistical laws 
emerge. These macroscopic laws arc the outcomes of 
self-organization and interactions among constituents 
through, which cannot be explained by the sum of the 
microscopic behaviors of individuals. Statistical laws can 
be extracted from time series, which are the most usual 
records of observable quantities in real world. The fractal 
and multifractal nature of time series have been exten- 
sively studies for different systems 

For a nonstationary time series, the detrended fluc- 
tuation analysis (DFA) can be adopted to explore its 
long-range autocorrelations and multifractal fea- 

tures (7H9[ . Alternatively, the detrending moving average 
(DMA) method can also be used for fractal analysis [10 - 
[l^ or multifractal analysis (iBj . Numerical experiments 
on monofractal time series unveil that the performance 
of the DMA method is comparable to the DFA method 
with slightly different priorities under different situations 



|16Hl8l |. However, for multifractal time series, the multi- 
fractal detrending moving average (MFDMA) performs 
better than the multifractal detrended fluctuation anal- 
ysis (MFDFA) [ll] . In addition, we note that both the 
DFA and DMA algorithms can be extended from one di- 
mension to higher dimensions for fractal and multifractal 
analysis [ll[ii-[2l|. 

A complex system usually contains several observable 
variables that exhibit long-range dependence or multi- 
fractal nature. In turbulent flows, the velocity, tempera- 
ture and concentration fields are embedded in the same 
space as joint multifractal measures [23 - l26| , in which the 
scaling behavior of the joint moments of two joint multi- 
fractal measures fii and fi2 are investigated 



(1) 
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where s is the box size. This framework has also been 
applied to study the joint multifractal nature between 
topographic indices and crop yield in agronomy 2T\ . [28| , 
trading volume and volatility in stock markets [29 1, nitro- 
gen dioxide and ground- level ozone [30j . heart rate vari- 
ability and brain activity of healthy humans 13 1[. and 
wind patterns and land surface air temperature [32|- 
For two stationary time series {x{i)} and {y{i)} of the 
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same length, the time-lagged cross-correlation or covari- 
ance provides another example [33l - [36j . 

(2) 



Cis) = {x{t)y{t + s)). 



E,y{s) = ( X{t)-X{t) Y{t)-Y{t) 



For two nonstationary time scries {a;(i)} and {y(j)} of the 
same length, one can study the followin g cro ss-correlation 
function between two detrended series [37|: 

C^y{s) = ( \x{t) - X{t)\ \Y{t + s)- Y{t + s)] ) , (3) 

where X{t) and Y{t) are certain trend functions of X{t) 
and Y{t), respectively. The detrended cross-correlation 
analysis (DCCA) was introduced to investigate the long- 
range power-law cross-correlations between two nonsta- 
tionary time series [H, [s^ : 



(4) 



where X{t) and Y{t) are certain trend functions of X{t) 
and Y{t) specific to moving windows of size s, respec- 
tively. The DFA method is a special case of this DCCA 
method when X{t) = Y{t). The DCCA method studies 
the temporal (not the cross-sectional) properties of two 
nonstationary time series, which is similar to the instant 
cross-correlations [13, . The significance of the cross- 
correlation can be assessed by statistical tests [12, l43| . 
The DCCA method has been applied to study volume 
change and price change of the Standard and Poor's 
(S&P) 500 Index [il], volatilities of the Brazilian agrar- 
ian commodities and stocks [i^, traffic flows [i^, and 
self-affine time series of taxi accidents [13] . 

More generally, the multifractal detrneded cross- 
correlation analysis was introduced to investigate the 
multifractal nature in the long-range power-law cross- 
correlations between two nonstationary time series (48| . 
which recovers the MFDFA method when X{t) = Y{t). 
We call this method the MF-X-DFA for the reason that 
win be clear in Sec. HI Note that the MF-X-DFA method 
is relevant to the multifractal height cross-correlation 
analysis with differences Hi]. The MF-X-DFA method 
has been applied to temporal and spatial seismic data 
[soj . sunspot numbers and river flow fluctuations 51 



stock index returns [52|, [53[, price- volume relationships 
in agricultural commodity futures markets [s^ ] , and spot 
and futures markets of WTI crude oil [55 1 . 

In this work, we introduce a variant of the MF-X-DFA 
algorithm, termed multifractal detrending moving aver- 
age cross-correlation analysis (MF-X-DM A) , which com- 
bines the ideas of MFDMA and DCCA. The main dif- 
ference between MF-X-DFA and MF-X-DMA is that the 
latter adopts local moving average as the trend function. 
Since the MFDMA algorithm outperforms the MFDFA 
algorithm for multifractal time series, we expect that 
the MF-X-DMA algorithm will show advantages over the 
MF-X-DFA algorithm. Our numerical experiments and 
real-work data analysis confirm this conjecture. 

The paper is organized as follows. The paper is orga- 
nized as follows. Section In] describes a unified framework 



of the MF-X-DFA and MF-X-DMA algorithms. Sec- 
tion HIT] performs extensive numerical experiments using 
fractal and multifractal time series with known analyt- 
ical expressions (bivariate fractional Brownian motions, 
two-component ARFIMA processes and binomial mea- 
sures) to investigate the performance of the algorithms. 
In Sec. IIVI we apply the algorithms to daily stock index 
returns and volatilities. We discuss and summarize our 
findings in Sec. |V] 



II. MF-X-DMA AND MF-X-DFA 

Consider two stationary time series {x{i)} and {y{i)} 
of the same length M, where i — 1,2, - ■ ■ ,M. Without 
loss of generality, we can assume that these two time 
series have zero means. Each time series is covered with 
Ms = [M/s] non-overlapping boxes of size s. The profile 
within the vth box [Z„ -I- 1 , -I- s] , where ly = (w — l)s, 
arc determined to be 

fe fc 
Xy{k) = xilv + j) and Yy{k) = y(^- + ■?■)' (5) 

where k ~ 1, ■ ■ ■ , s. Assume that the local trending func- 
tions of {Xy{k)} and {F„(fc)} are {Xy{k)} and {Yy{k)}, 
respectively. The cross-correlation for each box is calcu- 
lated as follows 

Fy{s) = -Y[Xvik) ~ Xy{k)\ (fc) - % (fc)] (6) 
The qth order cross-correlation is calculated as follows 

1/9 



Fxyiq,s) 

when g 7^ and 

Fxy{0,s) exp 



q/2 



v=l 



(7) 



2m ^-^ 



v=l 



We then expect the following scaling relation 



(8) 



(9) 



According to the standard multifractal formalism, the 
multifractal scaling exponent T(q) can be used to char- 
acterize the multifractal nature, which reads 



(10) 



where Df is the fractal dimension of the geometric sup- 
port of the multifractal measure Q . For time series anal- 
ysis, we have Df ^ 1. If the scaling exponent function 
T{q) is a nonlinear function of q, the signal has multifrac- 
tal nature. It is easy to obtain the singularity strength 
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function a (5) and the multifractal spectrum f{a) via the 
Legendre transform [56| 

( a^^yiq) =dT^y{q)/dq 

1 Uviq) = qotxy - Txy{q) 

There are many different methods for the determi- 
nation of Xy and Yy. The local detrending functions 
could be polynomials 0, Q , which recovers the MF-DXA 
method [48[. The local detrending function could also 
be the moving averages [lol . , in which the algorithm 
is called MF-X-DMA. To be more clear, we rename the 
MF-DXA algorithm as MF-X-DFA, and all muhifrac- 
tal analysis algorithms for cross-correlations based on 
local detrending are termed multifractal detrending/de- 
trended cross-correlation analysis (or MF-DCCA as ab- 
breviation). MF-X-DFA is an MF-DCCA method based 
on DFA and MF-X-DMA is an MF-DCCA method based 
on DMA. When X = Y, MF-X-DFA in Rcf. reduces 
to MF-DFA in Ref. [1], and MF-X-DMA reduces to MF- 
DMA in Ref. [l^. We note that the extension of the 
MF-DCCA algorithms to high dimensions is straightfor- 
ward d!]. 

The moving average function Z{t) oi Z G {X,Y} in a 
moving window [31 can be calculated as follows, 

-, r(n-i)(i-e)l 

J2 ^(^-^)' (12) 

where n is the window size, [gj is the largest integer 
smaller than g, \g~\ is the smallest integer larger than g, 
and 9 is the position parameter with the value varying 
in the range [0, 1]. Hence, the moving average function 
considers \{n — 1)(1 — 9)1 data points in the past and 
[{n — 1)9\ points in the future. We consider three spe- 
cial cases in this paper. The first case 9 = refers to 
the backward moving average [iTj . in which the moving 
average function Z{t) is calculated over all the past n—1 
data points of the signal. The second case 6* = 0.5 corre- 
sponds to the centered moving average [l^l, where Z{t) 
contains half past and half future information in each 
window. The third case = 1 is called the forward mov- 
ing average, where Z{t) considers the trend of n — 1 data 
points in the future. 



III. NUMERICAL EXPERIMENTS 

In order to investigate the validity and performance of 
the proposed MF-X-DCCA algorithms, we perform ex- 
tensive numerical experiments using bivariate fractional 
Brownian motions (bFBMs) j33l43a |. two-component 
autoregressive fractionally integrated moving average 
(ARFIMA) processes [H, [s^l, and binomial measures 
generated from the multiplicative p- model [ssj . By def- 
inition, there is no multifractality in bFBMs and two- 
component ARFIMA processes. Therefore, the hxy{q) 



function is independent of q and the r^y (q) function is lin- 
ear. In contrast, binomial measures is expected to possess 
multifractal nature. These three classes of time series are 
adopted to test the performance of the algorithms since 
the theoretical expressions of Hxx{q) are known for in- 
dividual time scries and we know the theoretical expres- 
sions of Hxy (q) for the first two classes and the numerical 
expression for the third class. Note that we have used H 
for theoretical values and h for estimated values below. 



A. Bivariate fractional Brownian motions 

A bivariate fractional Brownian motion [x{t), y{t)] with 
parameters {H^x, Hyy} G (0,1)^ is a self-similar Gaus- 
sian process with stationary increments, where x{t) and 
y{t) are two univariate fractional Brownian motions with 
Hurst indices H^x and Hyy and are the two components 
of the bFBM [ssl - fssj . The basic properties of multivari- 
ate fractional Brownian motions have been extensively 
studied [33l - l35| . Especially, it has been proven that the 
Hurst index H^y of the cross-correlation between the two 
components is [33l - l35l . |37| : 

Hxy = {HxX + Hyy)/2. (13) 

This property allows us to investigate the performances 
of the proposed algorithms on a solid foundation. 

An efficient simulation technique for univariate FBMs 
relies upon the embedding of the covariance matrix into 
a circulant matrix, whose square root can be easily 
obtained by the discrete Fourier transform [s^. This 
method is an exact simulation algorithm provided that 
the circulant matrix is definite positive. This algorithm 
can be generalized to simulate bivariate FBMs, which 
embeds the circulant of a block Toeplitz covariance ma- 
trix and use the fast Fourier transform to diagonalize the 
block circulant matrix [g^I • A detailed description of the 
simulation procedure can be found in Ref. [SJ, • 

In the simulation algorithm, the two Hurst indexes H^x 
and Hyy of the two univariate FBMs and their cross- 
correlation coefficient p are input arguments. We have 
generated a huge number of bFBMs, where H^x^ Hyy 
and p all vary from 0.1 to 0.9 with a spacing of 0.1. For 
a given triple of {Hxxi Hyy, p), 100 repeated simulations 
are conducted and 100 bFBMs with length of 2^^are gen- 
erated. In most cases, the definite positivity condition is 
not fulfilled. We then perform MF-X-DMA and MF-DFA 
on each bFBM to obtain the scaling exponent hxy The 
average over 100 repeated simulations is calculated. We 
have observed for each bFBM and each algorithm that 

hxy = {hxx + hyy)/2. (14) 

Our main findings are the following: (1) The exponent 
hxy is independent of the cross-correlation coefficient p; 
(2) The hxy{q) functions are independent of g, indicat- 
ing that the bFBMs are monofractals; (3) All the four 
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FIG. 1. (Color online) Multifractal detrended cross-correlation analysis of bivariate fractional Brownian motions. Comparisons 
are performed among three MF-X-DMA algorithms with 9 = 0, 0.5 and 1 and the MF-X-DFA algorithm. The results in plots 
(c-f) are averaged over 100 repeated simulations, (a) A typical example of bFBM with H^^ = 0.1, Hyy = 0.5 and p = 0.3. 
(b) Power-law dependence of the fluctuation functions Fxy{q,s) of the bFBM shown in plot (a) with respect to the scale s 
for q = —4, q = 0, and q = 4. The straight lines are the best power-law fits to the data. The results have been translated 
vertically for better visibility, (c) Scaling exponents hxy{q) with the theoretical values as a dashed line for Hxx = Hyy = 0.8 
and p = 0.5 (up panel) and H^x = 0.1, Hyy = 0.5 and p = 0.3 (low panel), (d) Independence of the scaling exponents h^y 
with respect to the cross-correlation coefficient p for Hxx = Hyy = 0.8 (up panel) and Hxx =0.1 and Hyy = 0.5 (low panel), 
(e) Differences Ahxy{q) between the estimated scaling exponents hxy and the theoretical exponents Hxy for bFBMs, where 
Hxx = 0.1 is fixed, Hyy varies from 0.1 to 0.9, and p takes different values, (e) Differences Ahxy{q) between hxy and Hxy for 
bFBMs with Hxx = Hyy varying from 0.1 to 0.9 and different p values. 



algorithms give nice estimates h^y of the scaling expo- 
nents, which are very close to the corresponding theo- 
retical Hxy values; and (4) The centered MF-X-DMA 
algorithm (9 = 0.5) and the MF-X-DFA algorithms have 
comparative performance and perform better than the 
backward and forward MF-X-DMA algorithms {9 = 
and 9 = 1). Since there are too many results to present 
in a concise way, we present a part of the results to man- 
ifest these findings. 

A typical example of the bFBM with Hxx = 0.1, 
Hyy = 0.5 and p = 0.3 is illustrated in Fig. [TJa) and 
the corresponding power-law dependence of the fluctua- 
tion functions Fxy{q, s) with respect to the scale s for the 
four algorithms is shown in Fig. [TJb). For MF-X-DMA 
algorithms, s should not be too large due to the finite- 
size effect. The scaling ranges span over two orders of 
magnitude for the MF-X-DMA algorithms and three or- 
ders of magnitude for the MF-X-DFA algorithm. In the 
determination of the scaling exponents h^y, we have used 
the same scaling ranges as in Fig.[ljb) for all the bFBMs 
and nice power-law relationships are observed. 

The two panels of Fig. [Ijc) show the hxy (q) functions 
for Hxx = Hyy = 0.8 and p = 0.5 and for Hxx = 0.1, 
Hyy = 0.5 and p = 0.3, respectively. Although there 



is a decreasing trend in each function, the theoretical 
functions Hxy{q) = 0.8 and Hxy{q) = 0.3 are well within 
the error bars, indicating that the hxy{q) functions are 
independent of the order q. Hence, the four algorithms 
arc able to correctly capture the monofractal nature of 
the bFBMs. We focus on g = 2 below. 

The two panels of Fig. [ijd) show the dependence of 
the scaling exponents hxy with respect to the cross- 
correlation coefficient p for Hxx = Hyy = 0.8 and for 
Hxx =0.1 and Hyy = 0.5, respectively. We find that 
the hxy functions (/i^=o, hl=°-^, hl=^ and h^^^) for the 
four algorithms are independent of p. This finding is 
very important since it distinguishes the temporal cross- 
correlations quantified by MF-DCCA algorithms and the 
cross-sectional correlation quantified by p. Two uncorre- 
lated time series may exhibit long-term power-law cross- 
correlation. In addition, the centered MF-X-DMA and 
the MF-X-DFA algorithms give similarly very accurate 
estimates of the scaling exponents with hxy^'^ ~ h^y^i 
Hxy = 0.8 for the up panel and Hxy = 0.3 for the low 
panel. In contrast, Fig.[TJc) shows that the backward and 
forward MF-X-DMA algorithms perform slightly worse 
and /if=o « hl=\ 

In order to compare the performance of the four algo- 
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rithms, we calculate the difference between the estimated 
exponent h^y and the theoretical exponent H^y: 



^hxy — hxy Hxy' 



(15) 



Figure [TJe) shows the dependence of Ahxy with respect 
to Hyy with a fixed H^x = 0.1 for p = 0.5 (top panel), 
p = 0.3 (middle panel) and p = 0.1 (bottom panel), while 
Fig. [IJf ) shows the dependence of A/i^jj, with respect to 
Hxx — Hyy for p = 0.9 (top panel), p = 0.5 (middle 
panel) and p — 0.1 (bottom panel). All the Ahxy val- 
ues in the two plots are less than 0.01, implying that all 
the four algorithms give good estimates. It is interest- 
ing to observe that /i^^" « h^y^ for all the cases. In 
addition, the centered MF-X-DMA and the MF-X-DFA 
algorithms outperform the backward and forward MF- 
X-DMA algorithms. We note that these conclusions also 
hold for other bFBMs. The relative performance between 
the centered MF-X-DMA and the MF-X-DFA algorithms 
are a little bit complicated. When Hj-x ^ Hyy, as shown 
in Fig.[IJe), the two algorithms have comparable perfor- 
mance since A/i^=°-5 « « 0. When H^^ = Hyy, 
as shown in Fig. ^i), the centered MF-X-DMA algo- 
rithm slightly outperforms the MF-X-DFA algorithm. In 
summary, the centered MF-X-DMA algorithm {9 — 0.5) 
is recommended for analyzing bivariate fractional Brow- 
nian motions. 



B. Two-component ARFIMA processes 

The power-law auto-correlations in stochastic variables 
can be modeled by an ARFIMA process (6l| : 



z{t) ^ Z{d,t)+e{t), 



(16) 



where d G (0, 0.5) is a memory parameter, is an in- 
dependent and identically distributed Gaussian variable, 
and 



Z[d, t) = an{d)z{t 



(17) 



in which a„((i) is the weight 

a„(d) = dr{n - d)/[T{l - d)T{n + 1)]. (18) 
The Hurst index Hzz is related to the memory parameter 

dby [lilel 



Hzz = 0.5 + d. 



(19) 



For the two-component ARFIMA processes discussed be- 
low, we take Z = X or Y. 

The two-component ARFIMA process is defined as fol- 
lows [131: 

r xit)=WX{d,,t) + {l~W)Y{d2,t) + ex{t) 

\y{t) = {l~W)X{di,t) + WYid2,t) + €y{t) ' ^^"^ 



where W £ [0.5,1] quantifies the coupling strength be- 
tween the two processes x{t) and y{t). When W = 1, 
x{t) and y{t) are fully decoupled and become two sep- 
arate ARFIMA processes as defined in Eq. ([TB)) . The 
cross-correlation between x jt) and y{t) increases when 
W decreases from 1 to 0.5 [53. To our knowledge, no 
general expression has been analytically derived for H^y 
When di > d2, the Hurst index H^x of x{t) decreases 
from 0.5 + di to certain value greater than 0.5 + d2 when 
W decreases from 1 to 0.5 [53|- In other words, H^x 
locates within the interval [0.5 -I- d2,0.5 + di]. When 
di = d2 = d, i.e. di d2, we obtain that 



Hyy 



0.5 



(21) 



which does not depend on the value of W . 

When IF = 1 and e^(t) = ty{t) = e{t), the two- 
component ARFIMA process becomes [s^ 



x{t) =X{di,t) + e{t) 
y{t)=Yid2,t) + eit) 



(22) 



If X and y are long-range power-law cross-correlated, it 
has been analytically derived that Eq. holds [4^ . 

The top panel (a-c) of Fig. [2] shows the results for the 
process in Eq. ^U^. Figure ^a) illustrates in log- log scale 
the dependence of the fluctuation functions Fxy{q, s) with 
respect to the scale s for 5 = —4, q = 0, and q = A 
for the process in Eq. pO| with di ~ d2 ~ 0.4. Nice 
power-law relations are observed, which are also evi- 
dent for other (^1,(^2) pairs. Figure [2Jb) shows the cor- 
responding scaling exponents hxy{q) for the process in 
Eq. pO| with di = d2 ^ 0.4. We note that the equation 



hxy{q) = [hxx{q) + hyy{q)]/2 holds for all the four curves. 
For the four algorithms, hxy{q) is close to the horizon- 
tal line H = 0.9, indicating that all the four algorithms 
correctly unveil the fractal nature of the two-component 
ARFIMA process. For g = 2, this plot shows that the 
MF-X-DFA gives the best estimate of hxy. Figure HJc) 
depicts the differences Ah^y between hxy and H^y with 
q = 2 for the process in Eq. (pil)) with different d values 
where di = d2 ~ d. It is found that, (1) the two MF-X- 
DMA algorithms with 9 = and 6 = 1 have the same 
performance, (2) the two MF-X-DMA algorithms with 
6* = and 6* = 1 perform better than the MF-X-DFA 
and the MF-X-DMA with 9 = 0.5 for relatively small 
d values, and (3) the two MF-X-DMA algorithms with 
9 = and 9 = 1 perform worse for large d values. 

The bottom panel (d-f) of Fig. [5] shows the results 
for the process in Eq. (|22p . Figure [H^d) illustrates in 
log-log scale the dependence of the fluctuation functions 
Fxy{q,s) with respect to the scale s for q = —4, q = 0, 
and q = 4: for the process in Eq. (|22p with di = 0.1 and 
d2 = 0.4. Nice power-law relations are observed, which 
are also evident for other (^1,^2) pairs. FigurejUJe) shows 
the corresponding scaling exponents h^yiq) for the pro- 
cess in Eq. with di =0.1 and (^2 = 0.4. Again, the 
equation hxy{q) = [hxx{q) + hyy{q)]/2 holds for all the 
four curves. For the four algorithms, hxy{q) is close to 



FIG. 2. (Color online) Multifractal detrended cross-correlation analysis of two-component ARFIMA processes. Comparisons are 
performed among three MF-X-DMA algorithms with 6 = 0, 0.5 and 1 and the MF-X-DFA method, (a) Power-law dependence 
of the fluctuation functions Fxy{q,s) with respect to the scale s for q = —4, q = 0, and q = 4 for the process in Eq. (|20|) 
with di = d2 = 0.4. The straight lines are the best power-law fits to the data. The results have been translated vertically for 
better visibility, (b) Scaling exponents hxy{q) for the process in Eq. (|20|l with d\ = d2 = 0.4. (c) Differences /^h^y between the 
estimated scaling exponents h^y and the theoretical exponents Hxy for the process in Eq. (|20p with different d values where 
d\ = d2 = d. id) Power-law dependence of the fluctuation functions Fxy{q,s) with respect to the scale s for q = —4, q = 0, 
and g = 4 for the process in Eq. (|22|) with di = 0.1 and d2 = 0.4. (e) Scaling exponents hxy{q) for the process in Eq. (|22p with 
di =0.1 and d2 = 0.4. (f) Differences Ah^y between h^y and H^y for the process in Eq. (|22p with different d values where 
di = d2 = d. 



the horizontal line H = 0.75, indicating that all the four 
algorithms correctly unveil the fractal nature of the two- 
component ARFIMA process. For q = 2, this plot shows 
that the MF-X-DMA algorithms with 6* = and 6* = 1 
give the best estimate of h^y Figure [^f) shows the dif- 
ferences between hj;y and H^y with q = 2 for the process 
in Eq. ((22)) with different d values where di = d2 = d. It 
is found that the MF-X-DFA and the MF-X-DMA with 
9 = 0.5 outperform the other two algorithms and give 
comparably nice estimates. 

He and Chen have investigated the two-component 
ARFIMA process defined in Eq. (PH)) of different lengths 
utihzing the DCCA method (the MF-X-DFA method 
with q = 2) and the DMCA method (the MF-X-DMA 
method with = and q = 2) and found that the DMCA 
method performs better in most cases and performs worse 
in a few cases However, both methods are prone to 
underestimate the exponents hxy f63| . Our results shown 
in Fig.^c) are consistent with their results for d = 0.15, 
0.25 and 0.35 in the sense that the MF-X-DMA method 
with 6* = outperforms the MF-X-DFA method. How- 
ever, we have obtained better estimates for hxy and there 
is no systematic underestimation. For instance, the three 
MF-X-DMA methods give h^y « 0.85 or Ah^y « when 
d = 0.35 as shown in Fig. [D^c). 

There are two subtle issues that might worsen the esti- 



mation of hxy As stated by Podobnik and Stanley [39j . 
they introduced a cutoff length M = 10* in their numeri- 
cal simulations and let the sum run from 1 to M, i.e., they 
set Oj = for j > Al. Our numerical experiments show 
that this cutoff seems optimal and a smaller or larger cut- 
off will worsen the estimation of the exponents. This find- 
ing applies both the MF-X-DMA and the MF-X-DFA. 
In addition, we stress that the upper bound of the scal- 
ing range for the MF-X-DMA algorithms should not be 
too large, because each moving averages are calculated 
within a window of size s. Let us take the MF-X-DMA 
algorithm with = as an example. In this case, the 
moving averages of the first s — 1 data points are not 
well defined. The bias becomes more significant for large 
window size s. When s is large, the Fxy{q,s) function 
bends downwards and the overall slope flattens. Similar 
arguments apply for other MF-X-DMA algorithms with 
different 6 values. In contrast, the MF-X-DFA algorithm 
does not suffer from this finite-size effect. 



C. Multifractal binomial measures 

We construct two binomial measures {x{i) : i = 
1,2,- ■• ,2'=} and {y{i) : i = 1, 2, • ■ • , 2'^'} from the p- 
model with known analytic multifractal properties as a 
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FIG. 3. (Color online) Multifractal detrended cross-correlation analysis of two cross-correlated binomial measures generated 
from the p-model with ~ 0.3 and Py = 0.4. Comparisons are performed among three MF-X-DMA algorithms with 6 = 0, 
0.5 and 1 and the MF-X-DFA method, (a) Power-law dependence of the fluctuation functions Fxy{q,s) with respect to the 
scale s for q = —4, q = 0, and q = 4. The straight lines are the best power-law fits to the data. The results have been 
translated vertically for better visibility, (b) Scaling exponents hxy{q) with the theoretical values as a dashed line. The insets 
show the hxy{q) curves and the corresponding [h^xiq) + hyy{q)]/2 curves, verifying the relation hxy{q) = [hxx{q) + hyy{q)\/2. 
(c) Multifractal mass exponents T(g) obtained from the MFDMA and MFDFA methods with the theoretical curve shown as a 
dashed line, (d) Multifractal spectra /(a) with respect to the singularity strength a for the four methods. The dashed curve 
is the theoretical multifractal spectrum, (e) Differences Ahxy{q) between the estimated mass exponents and their theoretical 
values for the three algorithms: MF-X-DFA, MF-X-DMA with e = and MF-X-DMA with 9 = 1. (f) Differences Ar(g) 
between the estimated mass exponents and their theoretical values for the three algorithms: MF-X-DFA, MF-X-DMA with 
e = and MF-X-DMA with 9 = 1. 



first example Each multifractal signal is obtained 

in an iterative way. We start with the zeroth itera- 
tion fc = 0, where the data set z(i) consists of one 
value, = 1. In the kth iteration, the data set 

{z^''^i) : i = 1, 2, • • • , 2'^} is obtained from 

for I = 1,2, ••• ,2*^^^. We notice that there are typos 
in the formula in Ref. [1^. When k oo, z^'''''{i) ap- 
proaches to a binomial measure, whose scaling exponent 
function Hzz{q) has an analytic form [sH fsl] 

HUq) = 1/9 - log2b^ + {l-p.Y]/q. (24) 
According to Eq. p^ . we have 

%M = ~\og^[vl + {l~VzY]- (25) 

In our simulation, we have performed fc = 16 iterations 
with px = 0.3 for x{i) and Py = 0.4 for y{i). The analytic 
scaling exponent functions H^xil) and Hyy{q) of x and 
y are expressed in Eq. (|24p . The two time series x and 
y are strongly correlated with a coefficient of 0.82, which 



is originated from the fact that the two sequences are 
constructed according to the same rules. The results are 
depicted in Fig. [31 

Figure ^b.) illustrates the power-law dependence of 
Fxy{q,s) against s for the four algorithms. Since the 
time series is not very long, we investigate —4 <q < 4 
to ensure the convergence of the qth moments [ij, [6^ . 
For the MF-X-DMA algorithms, there is a finite-size ef- 
fect since the moving averages at the ends of the time 
series are ill-defined. This effect becomes significant that 
deteriorates the estimation of F^y for large scales s. The 
scaling range is chosen as [2^,2"] for the three MF-X- 
DMA methods. In contrast, the MF-X-DFA method per- 
forms bad if the same scaling range is adopted. We use 
[2*^, 2^^] for the MF-X-DFA method which seems optimal. 
The algorithm-specific selection of the scaling range re- 
veals the difference in the applicability of the two types of 
methods. Figure EKa) shows that the power-law scaling 
is excellent for both positive and negative q values. 

The power-law scaling exponents {h^y^, h^y^, h^y^'^, 
and h^y^) are presented in Fig.[3^b), while the mass scal- 
ing exponents {t^^^, T^^y^, ^xy^'^ and rf~^) and the 
muhifractal spectra (/D/^^ fe=o^ f=i) are 
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illustrated in Fig. ^c) and Fig. [Sjd). It is evident that 
the MF-X-DMA method with 9 = 0.5 fails in a large 
part to correctly estimate the exponents, while the other 
three methods work much better. This finding is consis- 
tent with the conclusion that the MF-DMA method with 
6 = 0.5 performs much worse than the MF-DFA method 
and the MF-DMA methods with 6* = and 6* = 1 [H. 

The insets of Fig. [HJb) show an interesting feature for 
all the four algorithms that 

Ky{q) = [h:,Jq) + hyy{q)]/2, (26) 

no matter how accurate the estimates of an algorithm 
is. Similar relationship holds for individual monofractal 
ARFIMA signals [s^ and individual binomial measures 
[l5| . Hence, we can give the "theoretical" expression of 
Hxy{q) as follows 

Hxy{q)^[HxAq) + Hyyiq)]/2, (27) 

where H^xiq) and Hyy{q) are given in Eq. ([M)) . The 
theoretical line is plotted in Fig. ^h) as a dashed line. 
According to Eq. ([TU]), we obtain 

%yiq) = [%A<l)+Tyy{q)]/2, (28) 

which is shown in Fig.[3Ic) for comparison. Similarly, we 
have 

J-,y(a) = [Txx{a) + Tyy{a)\l2, (29) 

which is illustrated in Fig. EJd) as a dashed curve. 

In order to further assess the performance of the MF- 
X-DFA method and the two MF-X-DMA methods with 
= and = 1, we compare the empirical estimates of 
hxy{q) and Txy{q) with the theoretical values of Hxy{q) 
and Txyiq) by calculating the relative errors: 

^hxy{q) = hxy{q) - Hxy{q) (30) 

and 

^Txy{q) = Txyiq) - Txyiq)- (31) 

which are shown in Fig. ^e) and Fig. [3l[f). Roughly 
speaking, the MF-X-DMA algorithm with 6 = 1 performs 
best and the MF-DFA algorithm performs worst for most 
negative q values and the MF-X-DMA method with 9 = 
performs best and the MF-X-DMA method with 9 = \ 
performs worst for most negative q values. On average, 
the backward MF-X-DMA method (6* = 0) has the best 
performance and is thus recommended. 



IV. APPLICATION TO STOCK MARKET 
INDEX RETURNS AND VOLATILITIES 

We now apply the MF-X-DMA algorithms to investi- 
gate the temporal cross-correlations of the daily return 
and volatility time scries of the Dow Jones Industrial Av- 
erage (DJIA) and the National Association of Securities 



Dealers Automated Quotations (NASDAQ) index. The 
power-law cross-correlations between the DJIA volatil- 
ity and the NASDAQ volatility have been studied using 
the DCCA method ^ and the MF-X-DFA method [H. 
The closing prices of the DJIA and the NASDAQ from 
1971/2/5 to 2011/1/25 have been retrieved. The length 
of the time series is 10084. The return is defined as the 
daily difference of the logarithmic closing prices and the 
volatility is defined as the absolute value of the return. 

Figure mja) illustrates in log- log scale the dependence 
of the fluctuation functions Fxyiq, s) with respect to the 
scale s for q = —4, q = 0, and g = 4 for the returns. Ex- 
cellent power laws are observed spanning over two orders 
of magnitude for all the four algorithms. The resulting 
hxy functions are shown in Fig. IHb). All the four hxy 
functions arc monotonically decreasing, indicating that 
the cross-correlations between the index returns exhibit 
multifractal nature. We also find that 

Kfiq) > K=y\<l) > > (32) 

for —4 < (? < 4. When q = 2, all the hxy values are 
less than 0.6 and particularly « 0.5. This means 

that there is no significant linear long-term memory in 
the cross-correlations of the returns. The multifractal 
spectra fxyict) are plotted in Fig. Sfc). The singular- 
ity widths are all significantly greater than 0, confirming 
that the cross-correlations possess multifractal nature. 

Figure IH^d) illustrates in log-log scale the dependence 
of the fluctuation functions Fxyiq, ■^) with respect to the 
scale s for q = —4, q = 0, and q = 4 for the volatili- 
ties. Excellent power-law scaling is observed in the fluc- 
tuation functions for the MF-X-DMA algorithms with 
9 = and 1. However, for the MF-X-DFA algorithm 
and the MF-X-DMA algorithm with 9 = 0.5, there is 
a crossover in each curve. If we treat each curve with 
two scaling ranges and perform analysis on each scaling 
range, the resulting hxyiq) functions are not monotoni- 
cally decreasing and the multifractal spectra fxyia) are 
not concave. We thus focus on the two MF-X-DMA al- 
gorithms with 9 = and 1 that lead to one scaling range. 
The two hxyiq) functions are depicted in Fig. SJ^c). The 
two functions are monotonically decreasing and close to 
each other. When q = 2, hxy is close to 0.98, show- 
ing a very strong linear long-term memory in the cross- 
correlations between volatilities. We note that the rela- 
tion hxyiq) = [hxxiq) + hyyiq)]/2 does not hold, which 
is consistent with previous work [i^. Figure HJ^f) plots 
the two multifractal spectra. The large singularity width 
means that the cross-correlations between the two index 
volatilities exhibit multifractal nature. 

Our results for the volatility seems different from those 
in Refs. [H, IH. First of ah, the DJIA and NASDAQ 
time series are much longer in the current work. More 
importantly, the power-law scaling in the previous works 
exhibits significant fiuctuations j39l l48j , which makes it 
difficult to determine a proper scaling range. According 
to Fig.mjd), it is evident that the MF-X-DMA algorithms 
with 9 = and 1 significantly outperform the MF-X-DFA 
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FIG. 4. (Color online) Multifractal detrended cross-correlation analysis of the return time series (a-c) and the volatility time 
series (d-f ) for the DJIA index and the NASDAQ index. Comparisons are performed among three MF-X-DMA algorithms with 
^ = 0, 0.5 and 1 and the MF-X-DFA method. (a,d) Dependence of the fluctuation functions Fxy{q, s) with respect to the scale 
s for q — —4, q — 0, and q = 4. The straight lines are the best power-law fits to the data. The results have been translated 
vertically for better visibility. (b,e) Scaling exponents hxy{q) with respect to q. (c,f) Multifractal spectra ,fxy{a) with respect 
to the singularity strength a. 



algorithm. 



V. CONCLUSION AND DISCUSSION 

In this work, we have developed a class of MF-DCCA 
algorithms based on the detrending moving average anal- 
ysis. The performances of the MF-X-DMA algorithms 
are compared with the MF-X-DFA method by exten- 
sive numerical experiments on pairs of time series gen- 
erated from bivariatc fractional Brownian motions, two- 
component autorcgrcssive fractionally integrated moving 
average processes and binomial measures, which have 
theoretical expressions of the multifractal nature. In all 
cases, the scaling exponents h^y extracted from the MF- 
X-DMA and MF-X-DFA algorithms are very close to the 
theoretical values. 

For bivariate fractional Brownian motions, the scaling 
exponent hxy of the cross-correlation is found to be inde- 
pendent of the cross-correlation coefficient p between two 
time series. The MF-X-DFA and centered MF-X-DMA 
algorithms outperform the forward and backward MF- 
X-DMA algorithms. When H^^ ^ Hyy, the MF-X-DFA 
and centered MF-X-DMA algorithms show comparable 
performance. When H^x = Hyy, the centered MF-X- 
DMA algorithm performs slight better than the MF-X- 
DFA algorithm. Our numerical experiments verified the 
validity of the bFBM generating algorithm [13, IH, H^] . 
For two-component autoregressivc fractionally integrated 



moving average processes, we also found that the MF-X- 
DFA and centered MF-X-DMA algorithms have compar- 
ative performance, while the forward and backward MF- 
X-DMA algorithms perform slightly worse. All the four 
algorithms are able to correctly unveil the monofractal 
nature in the cross-correlations between the components 
of bFMBs and two-components ARFIMA processes. For 
binomial measures, the forward MF-X-DMA algorithm 
exhibits the best performance, the centered MF-X-DMA 
algorithms performs worst, and the backward MF-X- 
DMA algorithm outperforms the MF-X-DFA algorithm 
when the moment order g < and underpcrforms when 
q>0. 

In all the three mathematical models, the relation 
hxy = {hxx + hyy) 12 has been confirmed for all the four 
algorithms, where hxy, hxx and hyy are estimated scaling 
exponents. Previous works have shown that the MF-DFA 
and MF-DMA algorithms are able to give nice estimates 
for univariate signals, that is, hxx ~ Hxx and hyy « Hyy. 
It follows immediately that hxy ~ [Hxx + Hyy)/2. Com- 
bining the theoretical fact that Hxy ~ {Hxx + Hyy)/2, 
we obtain Hxy ~ hxy For monofractal time series, 
extensive numerical experiments unveiled that the per- 
formance of the DMA algorithms is comparable to the 
DFA algorithm and the centered DMA algorithm per- 
forms slightly better than DFA under certain situations 
[l6l - [T8| . This explains our numerical results on bFBMs 
and two-component ARFIMA processes. For multifrac- 
tal measures generated from the p-model, the backward 
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MF-DMA algorithm performs best [15|, which explains 
our findings for MF-DCCA algorithms. 

We applied these algorithms to the returns and volatil- 
ities of two US stock market indexes. For the returns, the 
centered MF-X-DMA algorithm gives the best estimates 
of hxy{q) since its hxy{2) is closest to 0.5, and the MF- 
X-DFA algorithm has the second best performance. For 
the volatilities, the forward and backward MF-X-DMA 
algorithms give similar results, while the centered MF- 
X-DMA and the MF-X-DFA algorithms fails to extract 
rational multifractal nature. These two applications are 
interesting since they showed that the choice of algo- 
rithms are automatic although we do not know which 
one should be used before analysis. The key message of 
our work is that we should use all the four algorithms 
and compare the results to make a choice. 
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Appendix A: Higher-dimensional MF-X-DMA 

In this work, we have focused on time series analysis. 
It is easy to generalize the one-dimensional MF-X-DMA 
algorithms to higher dimensions. The higher-dimensional 
MF-X-DMA algorithms are closely related to the MF-X- 
DFA algorithms [H], the DMA algorithms HO], and the 
MF-DMA algorithms [l^ in higher dimensions. 

Consider two physical quantities in d-dimension: 
{A(ii, • • • and • • • , irf)}, where ij = 

1,2, - ■■ ,Nj for j = 1,2,--- ,d. Before proceeding, we 
need to construct the difference matrixes {x{ii, ■ ■ ■ ,id)} 
and {y{ii, ■ • ■ ,id)} of X and Y. For simplicity, we de- 
note Z G {X, Y} and z as the corresponding difference 
matrix, which are related by the following equation: 



(Al) 



The matrix z is expressed as a square block matrix of 
size 2"^, whose block is z{ii, ■ ■ ■ ,id), where the intervals 
Ij = ij or [1 : ij — 1] for j = 1,2, •• • ,d. In real-world 
applications, we can focus on d ~ 2 and d ~ 3. 

For the two-dimensional case d = 2, the four blocks are 
12), z{ii - 1, 12), z{ii, 1 : Z2 - 1), and z(l : ii - 1, 1 : 



i2 — 1). According to Eq. (|A1[) . we have 

il-l 

Z{ii,i2) = z{ii,i2) + ^ z{j,i2) 

+x;^fe,j)+x!x!^(ji,j2). (A2) 

Since 

n — 1 12 — i 

^(■?'l'-'2) = - l.*2 - 1), 

Jl=l j2 = l 

11- l 

Y *2) = - 1, «2) - Z{ii - 1, 12 - 1), 

12- l 

z{i2,j) = Z{ii,i2 - 1) - Z{ii - l,i2 - 1), 

it follows that 

z{ii,i2) = Z{ii,i2) + Z{ii - 1,12 - 1) 

-Z(^l-l,^2)-Z(^l,^2-l), (A3) 

where Z{i,j) = if z x j = 0. 

For the three-dimensional case c? = 3, the eight blocks 
are z(ii, 12,23), z{l : ii - 1,12,23), z(«i,l : 22 - 1,^3), 
22, 1 : 23 - 1), z{l : 21 - 1, 1 : 22 - 1, 23), z{l : ii - 
1,22,1 : 23 - 1), 2(21,1 : 22 - 1,1 : 23 - 1), and z(l : 
2i — 1, 1 : 22 — 1, 1 : 23 — 1). We can derive similarly as 
the two-dimensional case that 

•z(«i,*2,23) = ^■(21,22,23) - Z{ii - 1,22 - 1,23 - 1) 

+Z{ii - 1, 22 - 1, 23) + ^■(21 - 1, 22, 23 - 1) 

-|-Z(2l,22 - 1,23 - 1) - Z(2l - 1,22,23) 
-Z(2l,22 - 1,23) - Z(2l,22,23 - 1), (A4) 

where Z{i,j,k) = if 2 x j x fc = 0. When 23 = 1, 
Eq. jXH reduces to Eq. ((A3|) . 

The algorithm of d-dimensional multifractal detrend- 
ing cross-correlation analysis is described as follows. 

Step 1. For each quantity z — x ov z — y, de- 
termine the moving averages Z{ii,--- ,id), where Sj ^ 
*j ^ ^ li'^j ^ l)^iJ ^i^d {^j} ^'''2 the position param- 
eters with the values varying in the range [0, 1] for j = 
1, 2, • • • , d. For each point located at (21, ■ • ■ , 2^) in the d- 
dimensional space, we extract a sub-matrix z(ki, ■ • ■ , kd) 
with size si x • • • x from the matrix z, where kj S 

- - 1)(1 - [{sj - 1)%J] ^ [mj,i,mj,2]- 

We calculate the cumulative sums Z'{ki, • • • , kd) of the 
points within the box: 

ki k2 

Z'(k,,--- ,kd)^ E •■■ E ziei,---Jd), (A5) 

and the moving average Z{ii,--- ,id) at location 
(21, ■ ■ ■ , id) is calculated as follows. 



mi,2 rnjo 



Z 



Si- ■ -Sd 



J2 ■■■ Z'{k,,.-- M)- (A6) 



fcl=nil,l fcl=md,l 
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Step 2. For each quantity, calculate the cumulative 
sums Q{ii,--- ,id) in a sliding window with size Si x 
• • • X Srf, where Sj ^ ij ^ Nj — [(sj — 1)6*^]. For each 
point located at (ii, • • • , id), we have 

il id 

Q= J2 E 4ki,---,kd). (A7) 

fcl=n-si+l kd=id-Sd + l 

Step 3. Detrend the matrix by removing the moving 
average function Ziii, ■ ■ ■ ,id) from Q{ii, ■ • ■ ,id), and ob- 
tain the residual matrix (11,12) as follows. 



e^(ii, • • • ,id) = Q{ii, ■■■ ,id) - Z{i 



ii • ■ • 



,u), (A8) 



where Sj ij ^ Nj - [(sj - l)9j\ . 

Step 4. Each residual matrix £2(11,- •• , irf) is parti- 
tioned into Ng-^ X ■ • ■ X Ns^ disjoint boxes of the same 
size si X ■ ■ ■ X Sd, where N^^ = [{Nj - Sj{l + 6j))/sj\. 
Each box can be denoted by ej^ ... for vj = 1, • • ■ , Ng. 
such that eJi^...^i,^(fci, • ■■ ,kd) = e'^^ly^ +fci, ■ • • ,lv^+kd) 
for 1 ^ fcj Sj, where ly- = VjSj. The cross-correlation 
between X and Y in each box is calculated as follows: 



, ■ ■ ■ 



1 



Kl = l Kd = l 



(fcl, • • • , fed) 



^4„-,Vd(ki,--- ,kd). (A9) 
S'fep 5. The gth order overall detrending cross- 



correlation function Fq{n) is calculated as follows: 



1 



E-Ei^^. 



19/2 



(AlO) 



vi — l 



Vd = l 



where = X]j=i ^j/*^ ^'^'^ 9 '^^^ take any real values 
except ioT q = 0. When g = 0, we have 



ln[Fo(n)] 



— ^ — F 

^ ^ Vi — 1 



^=d 

■E 

t'd=i 



ln|F,,^...,,J, (All) 



according to L'Hospital's rule. 

Step 6. Varying the box sizes Sj, we are able to de- 
termine the power-law relation between the fluctuation 
function Fq{s) and the scale s, 



F,(s) (A12) 



In real-world applications, one usually uses si = • • • = 
Sd = s. When Ng^ ^ {Nj — Sj{l + dj))/sj, one needs to 
start from different directions as for the DFA algorithm 
j66j or uses a random algorithm [63 ■ In addition, the 
box-by-box procedure is crucial for multifractal analy- 
sis, which was shown for high-dimensional MF-DFA [1^ 
and MF-DMA [Hi. However, the "traditional" proce- 
dure works well for high-dimensional fractals [20j . 



Appendix B: MATLAB codes for MF-X-DMA 



%% main function 

function [Fxxq, Fxyq , Fyyq , s] = ZQJIANG_MFXDMA_1D ( x , y, theta , q, s) 
if nargin < 5 

L = length (x) ; 
s [] ; i = 1.3; 
while round(10~i) L/4 
s = [ s round ( 1 " i ) ] ; 
i = i + 0.1; 

end 

clear L i 

end 

if nargin < 4, q = 2; end % DMA 
if nargin < 3, theta — . 5 ; end 
% X y must be transversal vector 
[i, j] = size(x);if i > 1 &&; j =^ 1 , x = x ' ; end 
[ i , j ] = size (y);if i > l&&:j = 1, y = y'; end 
clear i j 

Fxxq = zeros ( length ( q) , length(s)); 

Fxyq = Fxxq ; 

Fyyq = Fxxq; 

for i = l:length(s) 

x_re — myfun_MA ( X , s(i), theta); 
y_re = myfun_MA ( y , s(i), theta); 
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[Fxxq(:,i), Fxyq ( : , i ) , Fyyq ( : , i ) ] = myfun_Fq ( x_re , y.re , s(i), q) ; 

end 
end 

%% estimating the moving average residuals 

function ts_re = myfun_M A ( t s , s, theta) 

% estimate moveing average 

is = cumsum ( t s ) ; 

N = length ( ts ) ; 

A = zeros (s, N— s + 1); 

for k = 1 : s 

A(k, :) = ts (k:N-s+k) ; 

end 

MA = mean(A) ; 
clear A 

% moving average 

ts_re = ts(l + floor ((s-l)*thcta ): length ( ts )-ceil (( s-l)=i=(l- theta )) ) - IvlA; 

% If the residuals can not be completely covered by the series , we can cover the 

series from both sides. 
N= length ( ts_re ) ; 
n = fix(N/s) ; 
Is = N— n * s ; 
if Is ~= 
t s _ r c 1 
ts_re 

end 
end 

%% estimating the fluctuation scaling function 
function [Fxxq, Fxyq, Fyyq] = myfun_Fq ( x_re , y.re , s , q) 
n = length ( x_re )/ s ; 
X = reshape (x_re , s , n) ; 
Y = reshape(y_re , s , n) ; 
Fxx = mean( abs (X) . * abs (X) ) 
Fxy = mean( abs (X) . * abs (Y) ) 
Fyy = mean( abs (Y) . * abs (Y) ) 



ts. 
ts. 



rc { 
rc 1 



l:n*s); ts_re2 
ts_rc2 ] ; 



ts_re ( ls+l:N) ; 



Fxxq = zeros (length (q) 



1) 
1) 
1) 



Fxyq = zeros (length (q) 
Fyyq = zeros (length (q) 
for k = l:length(q) 

if q(k) ~= 

Fxxq ( k , 1 ) ~ (mean ( Fxx . 

Fxyq(k,l) = (mean(Fxy. 

Fyyq(k,l) = (mean(Fyy. 
elseif q(k) = 

Fxxq(k,l) 

Fxyq(k,l) 

Fyyq(k,l) 

end 

end 
end 



(q(k)/2))) 
(q(k)/2))) 
(q(k)/2))) 



(i/q(k)) 
(i/q(k)) 

(l/q(k)) 



exp (0.5 *mean( log ( Fxx) ) ) ; 
exp (0.5 *mean( log ( Fxy) ) ) ; 
exp (0.5 *mean( log ( Fyy) ) ) ; 
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Appendix C: MATLAB codes for MF-X-DFA 



%% main function 

function [Fxxq, Fxyq , Fyyq , s] = ZQ JIANG_MFXDFA_1D ( x , y, q, order, s) 
if nargin < 5 

L = length (x) ; 
s = [] ; i = 1.3; 
while round(10"i) <= L/4 
s = [ s round ( 1 " i ) ] ; 
i = i + 0.1; 

end 

clear L i 

end 

if nargin < 4, order = l;end 
if nargin < 3, q = 2; end % DFA 
x = cumsum(x) ; 
y = cumsum(y) ; 
N = length (x) ; 

Fxxq = zeros ( length ( q ) , length ( s ) ) ; 
Fxyq = zeros ( length ( q) , length ( s ) ) ; 
Fyyq — zeros ( length ( q) , length ( s ) ) ; 
for i = l:length(s) 

n = fix(N/s(i)); 

Is = N— n * s ( i ) ; 

if Is ~= 

xl = x(l:n*s(i)) ;x2 = x(ls+l:N); rsx = [xl x2 ] ; 
yl = y(l:n*s(i)) ;y2 = y(ls+l:N); rsy = [yl y2 ] ; 
n = 2 * n ; 

else 

rsx = X ; 
rsy = y; 

end 

X = reshape ( rsx , s ( i ) ,rL) ; 

Y = reshape ( rsy , s ( i ) ,rL) ; 
X = X' ; 

Y = Y' ; 

Fxx = zeros (1, n) ; 
Fxy = Fxx ; 
Fyy = Fxx; 
for j — 1 : n 

[Fxx(j) , Fxy(j) , Fyy(j)] = myfun.LocalResiduals ( 1 : s ( i ) , X( j , : ) , Y( j , : ) , 
order ) ; 

end 

for k = 1: length (q) 

if q(k) ~= 

Fxxq(k,i) = (mean(Fxx. ^ (q(k) /2) ) ) . ' (l/q(k) ) ; 

Fxyq(k,i) = (mean(Fxy."(q(k)/2))).'(l/q(k)); 

Fyyq(k,i) = (mean(Fyy . " ( q(k) /2) ) ) . ' ( 1 /q(k) ) ; 
elseif q(k) =^ 

Fxxq (k , i ) = exp (0.5 *mean( log (Fxx) ) ) ; 

Fxyq (k , i ) = exp (0.5 *mean( log (Fxy) ) ) ; 

Fyyq(k,i) = exp (0 . 5 *mean( log (Fyy ) ) ) ; 

end 

end 

end 
end 

%% Detrend the local trends 
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runction [ 


Kxx , 


Rxy , 


n 1 

R-yy] 


— 


r 

myiun. 


.LocalKcsi 


dual 


% n ^ 1 X 


= fx. 


'n; X 












% n = 2 X 


= [x. 


~ n; X 




V; 


X . ~ ( 


2)]; 




% n = 3 X 


= [x. 


' n ; X 


r(n- 


1); 


X . " ( n— 


2) ;x. ' (n- 


3)]; 


% n 4 X 


f 

= [x. 


' n ; X 


. (n- 


1); 


X . ( n— 


^J/z. ( n- 


3);x 


if n = 1 
















Z = [ t 


. * n ; t 


."(n 


-1)] ; 










elseif n= 


= 2 














Zj = [ t 


. n ; t 


• (n 


-Ij , t 


■ ( 


n-2) J ; 






elseif n= 


= 3 














Z = [t 


. * n ; t 


.'(n 


-l);t 




n-2);t 


."(n-3)]; 




elseif n = 


= 4 














Z = [t 


. ' n ; t 




-l);t 




n-2);t 


.^(ii-3);t 





-4)]; 



end 



A = x/Z; 
B = y/Z; 

Rxx = mean( abs (x— A*Z) . * abs (x— A^Z) ) 
Rxy = mean( abs (x— A=i=Z) . * abs (y— B^Z) ) 
Ryy = mean(abs (y— B*Z) . * abs (y— B*Z) ) 
end 
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